import numpy as np
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
import numpy as np
from scipy.interpolate import splrep, splev
#本题采用牛顿插值法实现
#定义差商函数
x = [2005,2010,2015,2020]
y = [72.8,74.2,75.2,76.4]
#y = [70.2,70.2,70.3,71.2]
#定义1~5阶差商函数
def diff1(i,j):
    return (y[i]-y[j])/(x[i]-x[j])
def diff2(i,j,k):
    return (diff1(i,j)-diff1(j,k))/(x[i]-x[k])
def diff3(i,j,k,l):
    return (diff2(i,j,k)-diff2(j,k,l))/(x[i]-x[l])
def diff4(i,j,k,l,m):
    return (diff3(i,j,k,l)-diff3(j,k,l,m))/(x[i]-x[m])
def diff5(i,j,k,l,m,n):
    return (diff4(i,j,k,l,m)-diff4(j,k,l,m,n))/(x[i]-x[n])
#定义牛顿插值函数
def newton3(xvalue):
    return y[0]+diff1(0,1)*(xvalue-x[0])+diff2(0,1,2)*(xvalue-x[0])*(xvalue-x[1])+diff3(0,1,2,3)*(xvalue-x[0])*(xvalue-x[1])*(xvalue-x[2])
def newton4(xvalue):
    return y[0]+diff1(0,1)*(xvalue-x[0])+diff2(0,1,2)*(xvalue-x[0])*(xvalue-x[1])+diff3(0,1,2,3)*(xvalue-x[0])*(xvalue-x[1])*(xvalue-x[2])+diff4(0,1,2,3,4)*(xvalue-x[0])*(xvalue-x[1])*(xvalue-x[2])*(xvalue-x[3])
print("地区一",newton3(2000),newton3(2013),newton3(2018),newton3(2025))
y = [70.2,70.2,70.3,71.2]
print("地区二",newton3(2000),newton3(2013),newton3(2018),newton3(2025))
#print("地区二",newton3(2005),newton3(2010),newton3(2015),newton3(2020))
x = [2000,2005,2010,2015,2020,2025]
y = [71.8,72.8,74.2,75.2,76.4]
y = [71.8,72.8,74.2,75.2,76.4,newton4(2025)]
rn = diff5(0,1,2,3,4,5)*(2025-x[0])*(2025-x[1])*(2025-x[2])*(2025-x[3])*(2025-x[4])
print("第二问地区一",newton4(2013),newton4(2018),newton4(2025))
print("第二问地区一误差",rn/newton4(2025))
x = [2000,2005,2010,2015,2020,2025]
y = [69.6,70.2,70.2,70.3,71.2]
y = [69.6,70.2,70.2,70.3,71.2,newton4(2025)]
rn = diff5(0,1,2,3,4,5)*(2025-x[0])*(2025-x[1])*(2025-x[2])*(2025-x[3])*(2025-x[4])
print("第二问地区二",newton4(2013),newton4(2018),newton4(2025))
print("第二问地区二误差",rn/newton4(2025))

#第三问，使用scipy进行三次样条插值
x = [2000,2005,2010,2015,2020]
y = [71.8,72.8,74.2,75.2,76.4]
f = splrep(x,y,k=3)
xnew = [2013,2018,2025]
ynew = splev(xnew,f)
print("第三问地区一",ynew[0],ynew[1],ynew[2])
x = [2000,2005,2010,2015,2020]
y = [69.6,70.2,70.2,70.3,71.2]
f = splrep(x,y,k=3)
xnew = [2013,2018,2025]
ynew = splev(xnew,f)
print("第三问地区二",ynew[0],ynew[1],ynew[2])

                                                                         